function idat = fsb_correct_intensity3(idat,gradientval)

% FSB - DEV: lets you manually correct changes in intensity over an image in y direction
% when occurring in a linear fashion
%
% EXAMPLE:
% idat = fsb_correct_intensity3(idat,10,10)
%
% INPUT:
% idat:     4D image
% gradientval: steepness of gradient for correction
%
% OUTPUT:
% idat: intensity corrected image
%
% CALLED BY:
% FSB.m
%
% NOTES
% 
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
%
% $ Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~
%$ Check if gradient value input provided
%~~~~~~~~~~~~~~~~~~~~~~~
if nargin<2
    gradientval = [3 3 3];
end

disp(['Correction value is x:' num2str(gradientval(1)) ' | y: ' num2str(gradientval(2)) ' | z: ' num2str(gradientval(3))]);
if gradientval == [0 0 0];
    disp('Correction not applied, returning...');
    return
end

%~~~~~~~~~~~~~~~~~~~~~~~
% Do actual intensity correction
%~~~~~~~~~~~~~~~~~~~~~~~

mean_idat(1) = mean(idat(:));

dim = size(idat);
for x=1:3
    corr_gradient = ones(1,dim(x));
    corr_gradient_tmp = 1:dim(x);
    corr_gradient_tmp = corr_gradient_tmp/dim(x)*gradientval(x);
    corr_gradient = corr_gradient_tmp+corr_gradient;
    if max(corr_gradient)>min(corr_gradient);
        corr_gradient = corr_gradient/(max(corr_gradient)-min(corr_gradient))*1;
    end
    figure(1995);
    subplot(3,1,x)
    plot(corr_gradient);

    switch x
        case 1
            corr_gradient = fliplr(corr_gradient);
            corr_gradient_tmp = repmat(corr_gradient',[1,dim(2),dim(3)]);
            corr_gradient_3d(:,:,:,x) = corr_gradient_tmp;
        case 2
            %corr_gradient = fliplr(corr_gradient);
            corr_gradient_tmp = repmat(corr_gradient,[dim(1),1,dim(3)]);
            corr_gradient_3d(:,:,:,x) = corr_gradient_tmp;
        case 3
            corr_gradient_tmp = repmat(corr_gradient_tmp',[1,dim(1),dim(2)]);
            corr_gradient_tmp = shiftdim(corr_gradient_tmp,1);
            corr_gradient_3d(:,:,:,x) = corr_gradient_tmp;
    end
end

corr_gradient_3d = mean(corr_gradient_3d,4);

idat = double(idat).*corr_gradient_3d;

mean_idat(2) = mean(idat(:));
idat = idat*(mean_idat(1)/mean_idat(2));

end
